#include <iostream>
#include <cstring>
#include <fstream>
#include <sstream>

#include <cv.h>
#include <cvaux.h>
#include <highgui.h>

using namespace std;
using namespace cv;

#define GLCM_DIS 7  //灰度共生矩阵的统计距离
#define GLCM_CLASS 64 //计算灰度共生矩阵的图像灰度值等级化
#define GLCM_ANGLE_HORIZATION 0  //水平
#define GLCM_ANGLE_VERTICAL   1	 //垂直
#define GLCM_ANGLE_DIGONAL    2  //对角


int calGLCM(IplImage* bWavelet,int angleDirection,double* featureVector,double featureCount)
{
	int i,j;
	int width,height;

	if(NULL == bWavelet)
		return 1;

	width = bWavelet->width;
	height = bWavelet->height;

	double * glcm = new double[GLCM_CLASS * GLCM_CLASS];
	int * histImage = new int[width * height];

	if(NULL == glcm || NULL == histImage)
		return 2;

	//灰度等级化---分GLCM_CLASS个等级
	uchar *data =(uchar*) bWavelet->imageData;
	for(i = 0;i < height;i++){
		for(j = 0;j < width;j++){
			histImage[i * width + j] = (int)(data[bWavelet->widthStep * i + j] * GLCM_CLASS / 256);
		}
	}

	//初始化共生矩阵
	for (i = 0;i < GLCM_CLASS;i++)
		for (j = 0;j < GLCM_CLASS;j++)
			glcm[i * GLCM_CLASS + j] = 0;

	//计算灰度共生矩阵
	int w,k,l;
	//水平方向
	
	if(angleDirection == GLCM_ANGLE_HORIZATION)
	{
		#pragma omp parallel for
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				l = histImage[i * width + j];
				if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width)
				{
					k = histImage[i * width + j + GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
				}
				if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width)
				{
					k = histImage[i * width + j - GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
				}
			}
		}
	}
	//垂直方向
	else if(angleDirection == GLCM_ANGLE_VERTICAL)
	{
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				l = histImage[i * width + j];
				if(i + GLCM_DIS >= 0 && i + GLCM_DIS < height) 
				{
					k = histImage[(i + GLCM_DIS) * width + j];
					glcm[l * GLCM_CLASS + k]++;
				}
				if(i - GLCM_DIS >= 0 && i - GLCM_DIS < height) 
				{
					k = histImage[(i - GLCM_DIS) * width + j];
					glcm[l * GLCM_CLASS + k]++;
				}
			}
		}
	}
	//对角方向
	else if(angleDirection == GLCM_ANGLE_DIGONAL)
	{
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				l = histImage[i * width + j];

				if(j + GLCM_DIS >= 0 && j + GLCM_DIS < width && i + GLCM_DIS >= 0 && i + GLCM_DIS < height)
				{
					k = histImage[(i + GLCM_DIS) * width + j + GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
				}
				if(j - GLCM_DIS >= 0 && j - GLCM_DIS < width && i - GLCM_DIS >= 0 && i - GLCM_DIS < height)
				{
					k = histImage[(i - GLCM_DIS) * width + j - GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
				}
			}
		}
	}
	double total=0;
	//normalization
	for (i = 0;i < GLCM_CLASS;i++)
	{
		for (j = 0;j < GLCM_CLASS;j++)
		{
			double cur = glcm[i * GLCM_CLASS + j];  	
			total += cur;
		}
	}
	for (i = 0;i < GLCM_CLASS;i++)
	{
		for (j = 0;j < GLCM_CLASS;j++)
		{
			glcm[i * GLCM_CLASS + j]=glcm[i * GLCM_CLASS + j]/total;  	
		}
	}
	//计算特征值
	double entropy = 0,energy = 0,contrast = 0,homogenity = 0;
	double mu1=0,mu2=0,sigma1=0,sigma2=0,cosum=0,correlation=0;
	double dis=0,IDM=0,CPro=0;
	#pragma omp parallel for
	for (i = 0;i < GLCM_CLASS;i++)
	{
		for (j = 0;j < GLCM_CLASS;j++)
		{
			//熵
			if(glcm[i * GLCM_CLASS + j] > 0)
				entropy -= glcm[i * GLCM_CLASS + j] * log10(double(glcm[i * GLCM_CLASS + j]));
			//能量
			energy += glcm[i * GLCM_CLASS + j] * glcm[i * GLCM_CLASS + j];
			//对比度
			contrast += (i - j) * (i - j) * glcm[i * GLCM_CLASS + j];
			//一致性
			homogenity += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];
			
			dis += glcm[i * GLCM_CLASS + j] * fabs(i-j);
						
			IDM += 1/(1+(i-j)*(i-j)) * glcm[i * GLCM_CLASS + j];
		}
	}
	for (i = 0;i < GLCM_CLASS;i++)
	{
	        for (j = 0;j < GLCM_CLASS;j++)
		{
			double cur = glcm[i * GLCM_CLASS + j];
			mu1 += i*cur;
			mu2 += j*cur;
		}
		
	}
	for (i = 0;i < GLCM_CLASS;i++)
	{
	        for (j = 0;j < GLCM_CLASS;j++)
		{
			double cur = glcm[i * GLCM_CLASS + j];
			sigma1 += (i-mu1)*(i-mu1)*cur;
			sigma2 += (j-mu2)*(j-mu2)*cur;
		}
		
	}
	double sigma=sqrt(sigma1*sigma2);
	for (i = 0;i < GLCM_CLASS;i++)
	{
	        for (j = 0;j < GLCM_CLASS;j++)
		{
			double cur = glcm[i * GLCM_CLASS + j];
			// printf ("%f\t",cur);
			correlation += (i-mu1)*(j-mu2)/sigma*cur;
		}
		// printf("\n");
		
	}
	for (i = 0;i < GLCM_CLASS;i++)
	{
	        for (j = 0;j < GLCM_CLASS;j++)
		{
			double cur = glcm[i * GLCM_CLASS + j];
			double pos = i+j-mu1-mu2;
			CPro += pow(pos,4)*cur;
		}
		
	}
	
	//返回特征值
	i = 0;
	featureVector[i++] = entropy;
	featureVector[i++] = energy;
	featureVector[i++] = contrast;
	featureVector[i++] = homogenity;
	
	/*cout<< "mu1:" << mu1 << endl ;
	cout<< "mu2:" << mu2 << endl;
	
	cout<< "sigma1:" << sigma1 << endl;
	cout<< "sigma2:" << sigma2 << endl;*/
	
	
	/*cout<< "entropy:" << entropy << endl ;
	cout<< "energy:" << energy << endl;
	cout << "contrast:" << contrast << endl;
	cout << "homogenity:" << homogenity << endl;
	cout << "correlation:" << correlation << endl;
	cout << "IDM:" << IDM << endl;
	cout << "dissimilarity:" << dis << endl;
	cout << "clusterPro:" << CPro << endl;*/
	cout << featureCount++ << ":" << entropy << " ";
	cout << featureCount++ << ":" << energy << " ";
	cout << featureCount++ << ":" << contrast << " ";
	cout << featureCount++ << ":" << homogenity << " ";
	cout << featureCount++ << ":" << correlation << " ";
	cout << featureCount++ << ":" << IDM << " ";
	cout << featureCount++ << ":" << dis << " ";
	cout << featureCount++ << ":" << log10(CPro) << " ";
	
	delete[] glcm;
	delete[] histImage;
	return 0;
}
static void read_csv(const string& filename, vector<Mat>& images, vector<float>& labels, char separator = ';') {
    std::ifstream file(filename.c_str(), ifstream::in);
    if (!file) {
        string error_message = "No valid input file was given, please check the given filename.";
        CV_Error(CV_StsBadArg, error_message);
    }
    string line, path, classlabel;
    while (getline(file, line)) {
        stringstream liness(line);
        getline(liness, path, separator);
        getline(liness, classlabel);
        if(!path.empty() && !classlabel.empty()) {
            images.push_back(imread(path,CV_LOAD_IMAGE_COLOR));
            labels.push_back(atof(classlabel.c_str()));
	    //printf("PUSH:%d\n",atoi(classlabel.c_str()));
        }
    }
}

int main(int argc, char **argv)
{
  vector<Mat> images;
  vector<float> labels;
  read_csv(string(argv[1]),images,labels);
  double featureCount=1;
  cout << labels[0] << " ";
  for(int i=0;i<images.size();i++)
  {
    Mat dest;
    cvtColor(images[i],dest,CV_BGR2GRAY);
    IplImage *img=new IplImage(dest);
    double fVec[5];
    calGLCM(img,GLCM_ANGLE_HORIZATION,fVec,featureCount);
    featureCount+=8;
  }
  cout << endl;
}








